{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we generate data from a mixture of exponentials with no censoring (this could be extended to censoring, but for our application we don't observe censoring). Furthermore, the data is not observed directly, but our instrument has binned the data into integer buckets. For example, the actual value is 4.3, but our instrument's resolution is too poor to measure the decimal point. So 4.3 goes into the \"4.x\" bucket. \n",
    "\n",
    "This model is easy to handle in *lifelines*. Instead of worrying about the binning, we can treat the system as an interval-censoring problem. That is, if an observation landed in the $i$th bin, then we know the _true_ data point occured somewhere between the $i$th and $i+1$th bin. This is precisely interval censoring. \n",
    "\n",
    "We can use *lifelines* custom model semantics to create a mixture model as well. The true model is:\n",
    "\n",
    "$$S(t) = p \\exp\\left(-\\frac{t}{\\lambda_1}\\right) + (1-p)\\exp\\left(-\\frac{t}{\\lambda_2}\\right)$$\n",
    "\n",
    "Therefore the cumulative hazard is:\n",
    "\n",
    "$$H(t) = -\\log(S(t)) = -\\log\\left(p \\exp\\left(-\\frac{t}{\\lambda_1}\\right) + (1-p)\\exp\\left(-\\frac{t}{\\lambda_2}\\right)\\right) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import exponential\n",
    "from lifelines.fitters import ParametricUnivariateFitter\n",
    "import autograd.numpy as np\n",
    "np.random.seed(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = [exponential(20) for _ in range(10000)] + [exponential(40) for _ in range(500)]\n",
    "counts_obs = np.bincount(T)\n",
    "T_obs = np.arange(np.amax(T))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 100.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwoAAAE/CAYAAADmEDQWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgd0lEQVR4nO3df7RdZX3n8fenNwQbgYhCVZJoAoa2QSu0KWg7VVtRY7HEaVFDtWKHTkor1VZdNrYWbawdpA7VjrQlLZlaLEZEh95qHGpF7VgFcxEqJpZFiEgSUZSfIggEvvPH3qmHs+/NPbk/cm/g/Vrrrpz9/Njn2ffutXM+Zz/POakqJEmSJKnXD830ACRJkiTNPgYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJjzhJNid53kyPYyYl+a9Jtie5O8lxe9n3M0l+Y7rGNh2S/F2SPxmw7eIklWTOdI+r5zkf9eekpP2PQUHSfiXJjUlO7Ct7TZLP7d6uqmOq6jPj7Gefv1jcx94NnFlVB1XV1TM9mEeT0ULLIOekJM02BgVJmgazIIA8Fdg8w2MY0yz4/UzI/jpuSZoIg4KkR5zeuw5Jjk8ykuSuJN9Kcm7b7F/bf+9op+c8O8kPJXlrkq8nuSXJ3yeZ37PfV7d1tyb5o77neXuSS5J8IMldwGva5/5CkjuS3JzkfUnm9uyvkvx2kuuTfDfJO5IcleTz7Xgv7m3fd4yjjjXJgUnuBoaAf09ywxj9fybJpiR3tv/+TF+To5J8sR3HPyZ5fNvvMe0x3toe16YkT2zr5ie5oD3WnUn+JMlQW/eaJP+W5M+T3Aq8o+3/9J4xHZ7k3iQ/0m6/JMk1bbvPJ/mJnrbHJflS+3v7EPCYPZwPQ0neneQ7SbYBJ/XVP+wuVfu3/ED7ePedp9OT3ARc3pZ/OMk329/fvyY5pi1fDbwSeHN7Xv1T/3O0f6P3JPlG+/OeJAe2dc9LsiPJG9u/681Jfr1nbL+YZEt73DuTvGms45akyTIoSHqkey/w3qo6BDgKuLgtf0777+Pa6TlfAF7T/vw8cCRwEPA+gCTLgL+keRH4ZGA+sKDvuVYClwCPA/4BeBD4PeAw4NnA84Hf7uvzIuCngGcBbwbWAa8CFgFPB04d47hGHWtV3VdVB7VtnllVR/V3bF/0fxz4C+AJwLnAx5M8oafZq4H/1h7rrrYtwGntsS9q+54B3NvW/V3b9mnAccALgd61DicA24AnAmuBj/Yd38uBz1bVLWnWVawHfrN9nvOB4fZF9lzgUuBC4PHAh4FfGeP3BPDfgZe0Y1oOnLKHtmN5LvDjNH8vgE8AS4EfAb5E8/emqta1j89pz6tfGmVff0jz9z4WeCZwPPDWnvon8YPz63TgvCSHtnUXAL9ZVQfTnB+XT+BYJGkgBgVJ+6NL23eZ70hyB80L+LE8ADwtyWFVdXdVXbGHtq8Ezq2qbVV1N/AWYFWa6SanAP9UVZ+rqvuBs4Dq6/+Fqrq0qh6qqnur6qqquqKqdlXVjTQvdp/b1+ecqrqrqjYDXwH+uX3+O2lejI61EHlPYx3PScD1VXVhO7YPAv8B9L6ovbCqvlJV3wP+CHh5e3fgAZoX7k+rqgfbY7yrvavwi8DvVtX3quoW4M+BVT37/EZV/a/2Oe8FLuqr/9W2DGA1cH5VXdk+z/uB+2heYD8LOAB4T1U9UFWXAJv2cLwvb9tur6rbgP8xwO+o39vb47oXoKrWV9V3q+o+4O3AM9Nz92kcrwTWVtUtVfVt4I+BX+upf6Ctf6CqNgJ3Az/aU7csySFVdXtVfWkCxyJJAzEoSNofvbSqHrf7h+679L1OB44G/qOdJvOSPbQ9Avh6z/bXgTk074AfAWzfXVFV9wC39vXf3ruR5OgkH2unqNwF/CnN3YVe3+p5fO8o2wcxuj2NdTz9fXf3771Dsr2v7gCasV8IXAZsaKfNnJPkAJo1EQcAN/cEuPNp3nEfbZ8AnwbmJTkhyWKad9j/T1v3VOCNfYFwUTv2I4CdVdUb1PqPp/94+49nb/1n/3Yq09lJbmj/rje2Vf1/2z2Np/9vd0TP9q1Vtatn+x5+cB78Ck0g+3qSzyZ59l4cgyTtFYOCpEe0qrq+qk6lecH6LuCSJI+lezcA4Bs0L1B3ewrNVJpvATcDC3dXJPlhmnfWH/Z0fdt/RfNO/dJ26tMfAJn40Qw81r3tu7v/zp7tRX11DwDfad/l/uOqWgb8DM2UnlfTvJC+DzisJ8QdUlXH9OznYb+fqnqQZirYqe3Px6rqu231duCdvYGwqua1dz9uBhYk6f1dPmUPx3vzKMfT63vAvJ7tJ42yj96x/yrNNLMTaaYILW7LM0rb0Yz2t/vGOH2aHVdtqqqVNOfzpfxgKp0kTTmDgqRHtCSvSnJ4VT0E3NEWPwR8u/33yJ7mHwR+L8mSJAfR3AH4UPvu7iXAL6VZBDyXZrrJeC/6DwbuAu5O8mPAb03RYY031vFsBI5O8qtJ5iR5BbAM+FhPm1clWZZkHs16gkuq6sEkP5/kGe00pLtoAsRDVXUz8M/A/0xySJrF1kcl6Z9q1e8i4BU003Eu6in/G+CM9m5Dkjw2yUlJDga+QBOKXpfkgCS/TDPPfywXt20XtnP91/TVX0MzbeuAJIOsYTiYJhTdShMw/rSv/ls8/Lzq90HgrWkWbx9GM43tA+M8J0nmJnllkvlV9QDN7/+h8fpJ0kQZFCQ90q0ANqf5JKD3Aqva9QP3AO8E/q2d2vIsmsWzF9J8ItLXgO8DvwPQriH4HWADzTvUdwO30LxgHMubaN59/i7NC98PTeFxjTnW8VTVrTR3At5I82L3zcBLquo7Pc0upFmc/E2aTxR6XVv+JJrQdBfwVeCzbVto7izMBbYAt7ftnjzOWK6keUf/CJo1GbvLR2gWIb+v3ddWmsXbtGtEfrndvo0maHx0D0/zNzTTpf6dZuFxf9s/olnofjvNeoGL2LO/p5kutJPmWPvXvVxAs47gjiSXjtL/T4AR4MvAte2YBvqyOJq1DDe2U57OoAlYkjQt8vApnpKkQbTv4t9BM63oazM8HEmSppx3FCRpQEl+Kcm8do3Du2neDb5xZkclSdL0MChI0uBW0iw6/QbNZ+ivKm/LSpIeoZx6JEmSJKnDOwqSJEmSOgwKkiRJkjrmzPQA+h122GG1ePHimR6GJEmStF+76qqrvlNVh0+0/6wLCosXL2ZkZGSmhyFJkiTt15J8fTL9nXokSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKnDoCBJkiSpw6AgSZIkqcOgIEmSJKljzkwPoN+1O+9k8ZqPj1l/49kn7cPRSJIkSY9O3lGQJEmS1DFQUEiyIsl1SbYmWbOHdr+SpJIs7yl7S9vvuiQvmopBS5IkSZpe4049SjIEnAe8ANgBbEoyXFVb+todDLweuLKnbBmwCjgGOAL4lyRHV9WDU3cIkiRJkqbaIHcUjge2VtW2qrof2ACsHKXdO4B3Ad/vKVsJbKiq+6rqa8DWdn+SJEmSZrFBgsICYHvP9o627D8l+UlgUVX1r0Iet2/bf3WSkSQjD95z50ADlyRJkjR9Jr2YOckPAecCb5zoPqpqXVUtr6rlQ/PmT3ZIkiRJkiZpkI9H3Qks6tle2JbtdjDwdOAzSQCeBAwnOXmAvpIkSZJmoUHuKGwCliZZkmQuzeLk4d2VVXVnVR1WVYurajFwBXByVY207VYlOTDJEmAp8MUpPwpJkiRJU2rcOwpVtSvJmcBlwBCwvqo2J1kLjFTV8B76bk5yMbAF2AW81k88kiRJkma/gb6Zuao2Ahv7ys4ao+3z+rbfCbxzguOTJEmSNAP8ZmZJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQMFhSQrklyXZGuSNaPUn5Hk2iTXJPlckmVt+eIk97bl1yT566k+AEmSJElTb854DZIMAecBLwB2AJuSDFfVlp5mF1XVX7ftTwbOBVa0dTdU1bFTOmpJkiRJ02qQOwrHA1uraltV3Q9sAFb2Nqiqu3o2HwvU1A1RkiRJ0r42SFBYAGzv2d7Rlj1MktcmuQE4B3hdT9WSJFcn+WySn5vUaCVJkiTtE1O2mLmqzquqo4DfB97aFt8MPKWqjgPeAFyU5JD+vklWJxlJMvLgPXdO1ZAkSZIkTdAgQWEnsKhne2FbNpYNwEsBquq+qrq1fXwVcANwdH+HqlpXVcuravnQvPkDDl2SJEnSdBkkKGwCliZZkmQusAoY7m2QZGnP5knA9W354e1iaJIcCSwFtk3FwCVJkiRNn3E/9aiqdiU5E7gMGALWV9XmJGuBkaoaBs5MciLwAHA7cFrb/TnA2iQPAA8BZ1TVbdNxIJIkSZKmzrhBAaCqNgIb+8rO6nn8+jH6fQT4yGQGKEmSJGnf85uZJUmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1GBQkSZIkdRgUJEmSJHUYFCRJkiR1DBQUkqxIcl2SrUnWjFJ/RpJrk1yT5HNJlvXUvaXtd12SF03l4CVJkiRNj3GDQpIh4DzgxcAy4NTeINC6qKqeUVXHAucA57Z9lwGrgGOAFcBftvuTJEmSNIsNckfheGBrVW2rqvuBDcDK3gZVdVfP5mOBah+vBDZU1X1V9TVga7s/SZIkSbPYnAHaLAC292zvAE7ob5TktcAbgLnAL/T0vaKv74IJjVSSJEnSPjNli5mr6ryqOgr4feCte9M3yeokI0lGHrznzqkakiRJkqQJGiQo7AQW9WwvbMvGsgF46d70rap1VbW8qpYPzZs/wJAkSZIkTadBgsImYGmSJUnm0ixOHu5tkGRpz+ZJwPXt42FgVZIDkywBlgJfnPywJUmSJE2ncdcoVNWuJGcClwFDwPqq2pxkLTBSVcPAmUlOBB4AbgdOa/tuTnIxsAXYBby2qh6cpmORJEmSNEVSVeO32ocOfPLSevJp7xmz/sazT9p3g5EkSZL2U0muqqrlE+3vNzNLkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeoYKCgkWZHkuiRbk6wZpf4NSbYk+XKSTyV5ak/dg0muaX+Gp3LwkiRJkqbHnPEaJBkCzgNeAOwANiUZrqotPc2uBpZX1T1Jfgs4B3hFW3dvVR07tcOWJEmSNJ0GuaNwPLC1qrZV1f3ABmBlb4Oq+nRV3dNuXgEsnNphSpIkSdqXBgkKC4DtPds72rKxnA58omf7MUlGklyR5KWjdUiyum0z8uA9dw4wJEmSJEnTadypR3sjyauA5cBze4qfWlU7kxwJXJ7k2qq6obdfVa0D1gEc+OSlNZVjkiRJkrT3BrmjsBNY1LO9sC17mCQnAn8InFxV9+0ur6qd7b/bgM8Ax01ivJIkSZL2gUGCwiZgaZIlSeYCq4CHfXpRkuOA82lCwi095YcmObB9fBjws0DvImhJkiRJs9C4U4+qaleSM4HLgCFgfVVtTrIWGKmqYeDPgIOADycBuKmqTgZ+HDg/yUM0oeTsvk9LkiRJkjQLDbRGoao2Ahv7ys7qeXziGP0+DzxjMgOUJEmStO/5zcySJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgwKkiRJkjoMCpIkSZI6DAqSJEmSOgYKCklWJLkuydYka0apf0OSLUm+nORTSZ7aU3dakuvbn9OmcvCSJEmSpse4QSHJEHAe8GJgGXBqkmV9za4GllfVTwCXAOe0fR8PvA04ATgeeFuSQ6du+JIkSZKmwyB3FI4HtlbVtqq6H9gArOxtUFWfrqp72s0rgIXt4xcBn6yq26rqduCTwIqpGbokSZKk6TJIUFgAbO/Z3tGWjeV04BMT7CtJkiRpFpgzlTtL8ipgOfDcvey3GlgNMHTI4VM5JEmSJEkTMMgdhZ3Aop7thW3ZwyQ5EfhD4OSqum9v+lbVuqpaXlXLh+bNH3TskiRJkqbJIEFhE7A0yZIkc4FVwHBvgyTHAefThIRbeqouA16Y5NB2EfML2zJJkiRJs9i4U4+qaleSM2le4A8B66tqc5K1wEhVDQN/BhwEfDgJwE1VdXJV3ZbkHTRhA2BtVd02LUciSZIkacoMtEahqjYCG/vKzup5fOIe+q4H1k90gJIkSZL2Pb+ZWZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUodBQZIkSVKHQUGSJElSx0BBIcmKJNcl2ZpkzSj1z0nypSS7kpzSV/dgkmvan+GpGrgkSZKk6TNnvAZJhoDzgBcAO4BNSYaraktPs5uA1wBvGmUX91bVsZMfqiRJkqR9ZdygABwPbK2qbQBJNgArgf8MClV1Y1v30DSMUZIkSdI+NsjUowXA9p7tHW3ZoB6TZCTJFUleujeDkyRJkjQzBrmjMFlPraqdSY4ELk9ybVXd0NsgyWpgNcDQIYfvgyFJkiRJ2pNB7ijsBBb1bC9sywZSVTvbf7cBnwGOG6XNuqpaXlXLh+bNH3TXkiRJkqbJIEFhE7A0yZIkc4FVwECfXpTk0CQHto8PA36WnrUNkiRJkmancYNCVe0CzgQuA74KXFxVm5OsTXIyQJKfTrIDeBlwfpLNbfcfB0aS/DvwaeDsvk9LkiRJkjQLDbRGoao2Ahv7ys7qebyJZkpSf7/PA8+Y5BglSZIk7WN+M7MkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpA6DgiRJkqQOg4IkSZKkDoOCJEmSpI45Mz2AvbV4zcfHrLvx7JP24UgkSZKkRy7vKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6hgoKCRZkeS6JFuTrBml/jlJvpRkV5JT+upOS3J9+3PaVA1ckiRJ0vQZNygkGQLOA14MLANOTbKsr9lNwGuAi/r6Ph54G3ACcDzwtiSHTn7YkiRJkqbTIHcUjge2VtW2qrof2ACs7G1QVTdW1ZeBh/r6vgj4ZFXdVlW3A58EVkzBuCVJkiRNo0GCwgJge8/2jrZsEAP1TbI6yUiSkQfvuXPAXUuSJEmaLrNiMXNVrauq5VW1fGje/JkejiRJkvSoN0hQ2Aks6tle2JYNYjJ9JUmSJM2QQYLCJmBpkiVJ5gKrgOEB938Z8MIkh7aLmF/YlkmSJEmaxcYNClW1CziT5gX+V4GLq2pzkrVJTgZI8tNJdgAvA85PsrntexvwDpqwsQlY25ZJkiRJmsXmDNKoqjYCG/vKzup5vIlmWtFofdcD6ycxRkmSJEn72KxYzCxJkiRpdjEoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6jAoSJIkSeowKEiSJEnqMChIkiRJ6hgoKCRZkeS6JFuTrBml/sAkH2rrr0yyuC1fnOTeJNe0P389xeOXJEmSNA3mjNcgyRBwHvACYAewKclwVW3paXY6cHtVPS3JKuBdwCvauhuq6tipHbYkSZKk6TTIHYXjga1Vta2q7gc2ACv72qwE3t8+vgR4fpJM3TAlSZIk7UuDBIUFwPae7R1t2ahtqmoXcCfwhLZuSZKrk3w2yc9NcrySJEmS9oFxpx5N0s3AU6rq1iQ/BVya5Jiququ3UZLVwGqAoUMOn+YhSZIkSRrPIHcUdgKLerYXtmWjtkkyB5gP3FpV91XVrQBVdRVwA3B0/xNU1bqqWl5Vy4fmzd/7o5AkSZI0pQYJCpuApUmWJJkLrAKG+9oMA6e1j08BLq+qSnJ4uxiaJEcCS4FtUzN0SZIkSdNl3KlHVbUryZnAZcAQsL6qNidZC4xU1TBwAXBhkq3AbTRhAuA5wNokDwAPAWdU1W3TcSCSJEmSps5AaxSqaiOwsa/srJ7H3wdeNkq/jwAfmeQYJUmSJO1jfjOzJEmSpA6DgiRJkqQOg4IkSZKkjun+HoV9avGaj0+o341nnzTh/Y7XV5IkSdofeUdBkiRJUodBQZIkSVKHQUGSJElSh0FBkiRJUscjajHzRE10EbQkSZL0SOUdBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdBgVJkiRJHQYFSZIkSR0GBUmSJEkdc2Z6APu7xWs+PmbdjWefNKF+45nofvfUbzzTtV9JkiTNTt5RkCRJktThHYVpNJm7BjOx3/3JdP0OvDsyOz1a7mjtT8e5P41VkjQx3lGQJEmS1DFQUEiyIsl1SbYmWTNK/YFJPtTWX5lkcU/dW9ry65K8aArHLkmSJGmajDv1KMkQcB7wAmAHsCnJcFVt6Wl2OnB7VT0tySrgXcArkiwDVgHHAEcA/5Lk6Kp6cKoPROMbb7rOTEwXmG3TqGZiPNO1yHyizzldC+33ZCam6c223/tknnO2TfWZiXNoPNPxwROPlr/n/sbf7eTMxO/vkXStnqjZeN2Ewe4oHA9sraptVXU/sAFY2ddmJfD+9vElwPOTpC3fUFX3VdXXgK3t/iRJkiTNYoMEhQXA9p7tHW3ZqG2qahdwJ/CEAftKkiRJmmVSVXtukJwCrKiq32i3fw04oarO7GnzlbbNjnb7BuAE4O3AFVX1gbb8AuATVXVJ33OsBla3m08HvjL5Q9Oj3GHAd2Z6ENqveQ5psjyHNFmeQ5qsH62qgyfaeZCPR90JLOrZXtiWjdZmR5I5wHzg1gH7UlXrgHUASUaqavmgByCNxvNIk+U5pMnyHNJkeQ5pspKMTKb/IFOPNgFLkyxJMpdmcfJwX5th4LT28SnA5dXcqhgGVrWfirQEWAp8cTIDliRJkjT9xr2jUFW7kpwJXAYMAeuranOStcBIVQ0DFwAXJtkK3EYTJmjbXQxsAXYBr/UTjyRJkqTZb6BvZq6qjcDGvrKzeh5/H3jZGH3fCbxzL8a0bi/aSmPxPNJkeQ5psjyHNFmeQ5qsSZ1D4y5mliRJkvToM9A3M0uSJEl6dJlVQSHJiiTXJdmaZM1Mj0ezX5JFST6dZEuSzUle35Y/Psknk1zf/nvoTI9Vs1uSoSRXJ/lYu70kyZXt9ehD7Yc5SKNK8rgklyT5jyRfTfJsr0PaG0l+r/1/7CtJPpjkMV6HNJ4k65Pc0n5Vwe6yUa89afxFez59OclPjrf/WRMUkgwB5wEvBpYBpyZZNrOj0n5gF/DGqloGPAt4bXverAE+VVVLgU+129KevB74as/2u4A/r6qnAbcDp8/IqLS/eC/wf6vqx4Bn0pxLXoc0kCQLgNcBy6vq6TQfHrMKr0Ma398BK/rKxrr2vJjmE0iX0nx/2V+Nt/NZExSA44GtVbWtqu4HNgArZ3hMmuWq6uaq+lL7+Ls0/zkvoDl33t82ez/w0hkZoPYLSRYCJwF/224H+AVg95dDeg5pTEnmA8+h+QRAqur+qroDr0PaO3OAH26/j2oecDNehzSOqvpXmk8c7TXWtWcl8PfVuAJ4XJIn72n/sykoLAC292zvaMukgSRZDBwHXAk8sapubqu+CTxxpsal/cJ7gDcDD7XbTwDuqKpd7bbXI+3JEuDbwP9up6/9bZLH4nVIA6qqncC7gZtoAsKdwFV4HdLEjHXt2evX2rMpKEgTluQg4CPA71bVXb117Zf/+fFeGlWSlwC3VNVVMz0W7bfmAD8J/FVVHQd8j75pRl6HtCftHPKVNKHzCOCxdKeTSHttstee2RQUdgKLerYXtmXSHiU5gCYk/ENVfbQt/tbu22ntv7fM1Pg06/0scHKSG2mmPP4CzXzzx7VTAMDrkfZsB7Cjqq5sty+hCQ5ehzSoE4GvVdW3q+oB4KM01yavQ5qIsa49e/1aezYFhU3A0naF/1yaRTzDMzwmzXLtXPILgK9W1bk9VcPAae3j04B/3Ndj0/6hqt5SVQurajHNdefyqnol8GnglLaZ55DGVFXfBLYn+dG26PnAFrwOaXA3Ac9KMq/9f233OeR1SBMx1rVnGHh1++lHzwLu7JmiNKpZ9YVrSX6RZq7wELC+/VZnaUxJ/gvw/4Br+cH88j+gWadwMfAU4OvAy6uqf7GP9DBJnge8qapekuRImjsMjweuBl5VVffN4PA0iyU5lmYx/FxgG/DrNG/GeR3SQJL8MfAKmk/zuxr4DZr5416HNKYkHwSeBxwGfAt4G3Apo1x72hD6PpppbfcAv15VI3vc/2wKCpIkSZJmh9k09UiSJEnSLGFQkCRJktRhUJAkSZLUYVCQJEmS1GFQkCRJktRhUJAkSZLUYVCQJEmS1GFQkCRJktTx/wHy/IQ07CNf5AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 936x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(13,5))\n",
    "plt.hist(counts_obs, bins=T_obs, density=True)\n",
    "plt.title(\"Histogram of observed durations\")\n",
    "plt.xlim(0, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.015428571428572\n"
     ]
    }
   ],
   "source": [
    "# To help the model avoid the \"non-identibility\" problem, we can set the *upper* bound of the first lambda to be \n",
    "# the average of the observed data, and the *lower* bound of the second lambda to the same value. Why? \n",
    "# We'd like to partition the postive reals into two halves, each containing one of the lambdas. The sample\n",
    "# mean of the data is v = p * lambda1 + (1-p) * lambda2, which has the property lambda1 < v < lambda2, therefore  \n",
    "# it will partition the space correctly. \n",
    "mean_obs = np.average(T_obs, weights=counts_obs)\n",
    "print(mean_obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>model</th>\n",
       "      <td>lifelines.MixtureExponential</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>number of observations</th>\n",
       "      <td>10500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>number of events observed</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>log-likelihood</th>\n",
       "      <td>-42207.84</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hypothesis</th>\n",
       "      <td>lambda1 != 10.0077, lambda2 != 40.0309, p != 0.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div><table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>coef lower 95%</th>\n",
       "      <th>coef upper 95%</th>\n",
       "      <th>z</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>lambda1</th>\n",
       "      <td>18.04</td>\n",
       "      <td>1.34</td>\n",
       "      <td>15.41</td>\n",
       "      <td>20.67</td>\n",
       "      <td>5.98</td>\n",
       "      <td>&lt;0.005</td>\n",
       "      <td>28.77</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>lambda2</th>\n",
       "      <td>29.89</td>\n",
       "      <td>5.33</td>\n",
       "      <td>19.44</td>\n",
       "      <td>40.34</td>\n",
       "      <td>-1.90</td>\n",
       "      <td>0.06</td>\n",
       "      <td>4.13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>p</th>\n",
       "      <td>0.79</td>\n",
       "      <td>0.18</td>\n",
       "      <td>0.44</td>\n",
       "      <td>1.14</td>\n",
       "      <td>1.63</td>\n",
       "      <td>0.10</td>\n",
       "      <td>3.28</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table><div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>AIC</th>\n",
       "      <td>84421.68</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrrrrrr}\n",
       "\\toprule\n",
       "{} &  coef &  se(coef) &  coef lower 95\\% &  coef upper 95\\% &     z &    p &  -log2(p) \\\\\n",
       "\\midrule\n",
       "lambda1 & 18.04 &      1.34 &           15.41 &           20.67 &  5.98 & 0.00 &     28.77 \\\\\n",
       "lambda2 & 29.89 &      5.33 &           19.44 &           40.34 & -1.90 & 0.06 &      4.13 \\\\\n",
       "p       &  0.79 &      0.18 &            0.44 &            1.14 &  1.63 & 0.10 &      3.28 \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.MixtureExponential:\"MixtureExponential_estimate\", fitted with 10500 total observations, 10500 interval-censored observations>\n",
       "   number of observations = 10500\n",
       "number of events observed = 0\n",
       "           log-likelihood = -42207.84\n",
       "               hypothesis = lambda1 != 10.0077, lambda2 != 40.0309, p != 0.5\n",
       "\n",
       "---\n",
       "          coef   se(coef)   coef lower 95%   coef upper 95%     z      p   -log2(p)\n",
       "lambda1  18.04       1.34            15.41            20.67  5.98 <0.005      28.77\n",
       "lambda2  29.89       5.33            19.44            40.34 -1.90   0.06       4.13\n",
       "p         0.79       0.18             0.44             1.14  1.63   0.10       3.28\n",
       "---\n",
       "AIC = 84421.68"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class MixtureExponential(ParametricUnivariateFitter):\n",
    "\n",
    "    _fitted_parameter_names = ['lambda1', 'lambda2', 'p']\n",
    "    _bounds = [(0, mean_obs), (mean_obs, None), (0, 1)]\n",
    "\n",
    "    def _cumulative_hazard(self, params, t):\n",
    "        l1, l2, p = params\n",
    "        return -np.log(p * np.exp(-t / l1) + (1-p) * np.exp(-t / l2))\n",
    "\n",
    "\n",
    "model = MixtureExponential()\n",
    "model.fit_interval_censoring(\n",
    "    lower_bound=T_obs, \n",
    "    upper_bound=T_obs + 1, \n",
    "    weights=counts_obs, \n",
    "    initial_point=np.array([mean_obs / 2, mean_obs * 2, 0.5])\n",
    ")\n",
    "\n",
    "model.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAEvCAYAAACdahL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/jElEQVR4nO3dd3hc5Z3//fetMhr13oslF9wblm3AQCihhB4gAUKyAZKwubJZkuyGLMlv2WzyPNmH3SVk0/aXkGxC2MJCcABTQwsBQosNtiTLNu62uixZbTT93M8fIxvbuEojzYz0eV2XL0szozO3j4/00bnb11hrERERkYmVFOsGiIiITEUKYBERkRhQAIuIiMSAAlhERCQGFMAiIiIxoAAWERGJgZSJfLOioiJbW1s7kW8pIiISM+vWrdtnrS0+2nMTGsC1tbWsXbt2It9SREQkZowxu4/1nLqgRUREYkABLCIiEgMKYBERkRiY0DFgEZEDgsEgLS0t+Hy+WDdFZMzcbjdVVVWkpqae9NcogEUkJlpaWsjOzqa2thZjTKybIzJq1lp6enpoaWmhrq7upL9OXdAiEhM+n4/CwkKFryQ8YwyFhYWn3JujABaRmFH4ymQxmmtZASwiIhIDCmARmbKMMXz6058++HkoFKK4uJgrrrgCgDVr1nDPPfcc9xgPPPAAbW1tUWvTAw88QHFxMUuWLDn4p7m5OWrHj6Zdu3bxP//zPwc/X7t2LXfcccdxv+aVV145eH7Hqq+vj3//938/+HlbWxvXX399VI79+OOPj/t5VwCLyJSVmZlJU1MTXq8XgBdeeIHKysqDz1911VXcddddxz3GaAI4FAod9/kbbriB9evXH/wzb968Uzr+RDkygOvr6/nRj340Ye9/ZABXVFTw6KOPRuXYCmARkXF22WWX8fTTTwPw0EMPcdNNNx187oEHHuDLX/4yAFdffTUPPvggAD//+c+5+eabefTRR1m7di0333wzS5Yswev1Ultby759+4DIHeF5550HwD/+4z/ymc98hlWrVvGZz3yG7u5urrvuOpYvX87y5cv505/+dNx2PvbYY1x44YVYa2lvb+e0006jo6ODBx54gKuvvprzzjuPWbNm8Z3vfOfg19x3330sWLCABQsW8G//9m9AJDTnzp3LF77wBebPn8/FF1988BeQ7du3c+mll7Js2TLOOeccNm/eDMAtt9zCHXfcwVlnncX06dMPhtxdd93Fa6+9xpIlS/jBD35w2N3tO++8w5lnnsnSpUs566yz2LJly0n9f3g8Hm677TZWrFjB0qVLeeKJJwDYuHEjK1asYMmSJSxatIitW7dy1113sX37dpYsWcKdd97Jrl27WLBgwcH/u2uuuYaLLrqI2tpafvKTn3DfffexdOlSzjjjDHp7ewH4xS9+wfLly1m8eDHXXnsdfQODvPHGG6xZs4Y777yTJUuWsH379mOem7HQMiQRibnvPLmR5raBqB5zXkUO375y/glfd+ONN/Ld736XK664goaGBm677TZee+21D73u/vvvZ9WqVdTV1fH973+ft956i4KCAn7yk59w7733Ul9ff8L3am5u5vXXXyc9PZ1PfepTfO1rX+Pss89mz549XHLJJWzatAmAhx9+mNdff/3g17355pt8/OMfZ/Xq1fz0pz/lueee4zvf+Q5lZWVAJOyamprIyMhg+fLlXH755Rhj+PWvf83bb7+NtZaVK1fykY98hPz8fLZu3cpDDz3EL37xCz75yU+yevVqPv3pT3P77bfzs5/9jFmzZvH222/zpS99iZdffhmA9vZ2Xn/9dTZv3sxVV13F9ddfzz333MO9997LU089BUS6lw+YM2cOr732GikpKbz44ot861vfYvXq1Sc8R9/73ve44IIL+NWvfkVfXx8rVqzgox/9KD/72c/4yle+ws0330wgECAcDnPPPffQ1NTE+vXrgcgvF4dqamrivffew+fzMXPmTP75n/+Z9957j6997Ws8+OCDfPWrX+Xaa6/l85//PP6Qw9///d9z/y//g2/8zVe56qqruOKKKw52aV944YXHPDejpQAWkSlt0aJF7Nq1i4ceeojLLrvsmK8rLS3lu9/9Lueffz6PPfYYBQUFp/xeV111Fenp6QC8+OKLh3VxDgwMMDQ0BES6oH/yk5986Ot//OMfs2DBAs4444zD7tQvuugiCgsLAbj22mt5/fXXMcbw8Y9/nMzMzIOPv/baa1x11VXU1dWxZMkSAJYtW8auXbsYGhrijTfe4BOf+MTB4/r9/oMfX3PNNSQlJTFv3jw6OztP+G/t7+/ns5/9LFu3bsUYQzAYPKlz9Pzzz7NmzRruvfdeILJcbc+ePZx55pl873vfo6WlhWuvvZZZs2ad8Fjnn38+2dnZZGdnk5uby5VXXgnAwoULaWhoAGD9hgbuvvtu+vr6GPZ4uOCjH/3QcU50bkZLASwiMXcyd6rj6aqrruLrX/86r7zyCj09Pcd8XWNjI4WFhccd801JScFxHIAPrQs9EIYAjuPw1ltv4Xa7T7qdLS0tJCUl0dnZieM4JCVFRhGPXAJzoiUxaWlpBz9OTk7G6/XiOA55eXkH7yaP9zXW2hO29e677z74y8quXbsOdsWfiLWW1atXM3v27MMenzt3LitXruTpp5/msssu4+c//znTp08/7rEObXNSUtLBz5OSkggGg3j8IW677Vb+5+FHWbx4Cf/1n7/h1Vf/+KHjnOjcjJbGgEVkyrvtttv49re/zcKFC4/5mnfeeYdnn32W9957j3vvvZedO3cCkJ2dzeDg4MHX1dbWsm7dOoDjdrlefPHF/PjHPz74+Yl+uIdCIW677TYeeugh5s6dy3333XfwuRdeeIHe3l68Xi+PP/44q1at4pxzzuHxxx9neHgYj8fDY489xjnnnHPM4+fk5FBXV8dvf/tbIBKEGzZsOG6bjvy3H6q/v//ghLYHHnjguMc51CWXXMKPf/zjgyH/3nvvAbBjxw6mT5/OHXfcwdVXX01DQ8Nx3/9YrLUEww6BsEMo7DA0NERFeQXBYJBH/veho/7bRnNuToYCWESmvKqqquMun/H7/XzhC1/gV7/6FRUVFXz/+9/ntttuw1rLLbfcwhe/+MWDk7C+/e1v85WvfIX6+nqSk5OPecwf/ehHrF27lkWLFjFv3jx+9rOfHXzu4YcfPmwZ0htvvME//dM/cc4553D22Wdz33338ctf/vLgmPGKFSu47rrrWLRoEddddx319fWcfvrp3HLLLaxYsYKVK1fy+c9/nqVLlx73PPz3f/83//Ef/8HixYuZP3/+wQlQx7Jo0SKSk5NZvHgxP/jBDw577hvf+Abf/OY3Wbp06QlnfR/q7rvvJhgMsmjRIubPn8/dd98NwCOPPMKCBQtYsmQJTU1N/MVf/AWFhYWsWrWKBQsWcOedd57w2KGww6AvRDDkYDAkJyVx9z/8I+efu4qLzj+XWad9cNd944038q//+q8sXbqU7du3n/K5ORnmZLoSoqW+vt6uXbt2wt5PROLXpk2bmDt3bqybkfAeeOAB1q5de9QxY4lwHIsvGMYfdkg2hqRjdNFbLI615Ka7RvU+R7umjTHrrLVHnaGnMWAREZmUrLUEQg7eYBiA1CQDxM/2pwpgEZEEdsstt3DLLbfEuhmn7Ne//jU//OEPD3ts1apV/PSnP43K8UNhB28gTMhaUpIMJo6C9wAFsIiITLhbb72VW2+9NerHdRyLLxTGH4p0N6cmxe9UJwWwiMSMtVYVkSQqDsxuHg7Eprt5NPOp4vdXAxGZ1NxuNz09PaP6wSVyqFDYYcgfwhMIk5xkSElKYqLDt6en55TWdIPugEUkRqqqqmhpaaG7uzvWTZEEFbnrtYQcB2PMmO8oLWAtpLuOvXzsWNxuN1VVVaf0NQpgEYmJ1NRU6urqYt0MSUCOY2nt89LY2kfYsRRkph1zadGpCDuWIX+QS+eWR6GVJ6YAFhGRhNE/HKShtY+eIT/5GWm4UhJ3JFUBLCIicc8fCrOta4htXUO4U5MpzUmPdZPGTAEsIiJxy1pLW5+XhtZ+wmFLUWYaSUmTY+a8AlhEROJSvzdIU2s/XYN+8jNSSUs/9clR8UwBLCIicSUQctjWNcTWzkHcqcmU5Zza8p5EoQAWEZG4YK2lo9/HhpY+gmFLUdbk6W4+GgWwiIjE3IAv0t3c2e8jP9NF7iTrbj4aBbCIiMRMIOSwvTvS3exKSaYsN/FnN5+sEy6gMsb8yhjTZYxpOuSxAmPMC8aYrSN/549vM0VEZDKJdDd7+cOWLrZ1DVGQmUZuemqsmzWhTmYF8wPApUc8dhfwkrV2FvDSyOciIiInNOgL8vbOXt7a0YMrOYmirDSSJ/FY77GcsAvaWvuqMab2iIevBs4b+fg3wCvA30WzYSIiMrkEww47uofY3DGIKzlpUmymMRajHQMutda2j3zcAZRGqT0iIjLJWGvpHPDR2NqPNximMHNq3vEeacyTsKy11hhzzHpixpjbgdsBampqxvp2IiKSQIb8IZpa+2nv85GXkUpx1tQa5z2e0e5i3WmMKQcY+bvrWC+01t5vra231tYXFxeP8u1ERCSRBMMO73cO8vKmTvZ7ApTmpOFOnfxLi07FaO+A1wCfBe4Z+fuJqLVIREQSlrWWrgEfDa39+AJhChKou3k4EKLXE5iw9zthABtjHiIy4arIGNMCfJtI8D5ijPkcsBv45Hg2UkRE4p/HH2Jj2wCt+73kpqdSlJ0Y3c29ngDPNbXz4qYuZpVk8amV0ybkfU9mFvRNx3jqwii3RUREElAo7LBzn4dN7QOkpiRRmpOGMfF/17u3d5inG9t5fds+HGtZUVvA+bMnbqhUO2GJiMioWGvpHvLTsLeP4UCY/EwXKUmjnVo0May1bGof4MmGdtbv7cOVnMSFc0q4bGE5RVlpDPmDE9YWBbCIiJwyjz9Ec3ukuznHnUpxnHc3hx3LOzt7ebqxje3dHnLcKXxiWRUXzSsl25168DUTSQEsIiInLRR22LXPQ3PHIClJhpLs+O5u9ofC/HFLN083ttM16Kcsx83nzq7j3FnFuFJie7euABYRkZPSNeijoaUPjz9MQZx3Nw/4gjy/sYPfb+xkyB9iVkkWn145jWXT8uOmxKECWEREjsvjD7GpfYCW/V6y3SmUZLtj3aRj6hzw8UxjO69s6SYQdji9Jp8rF5czuzQ77u7UFcAiInJUobDDrh4Pze3x3928vXuIpxraeHtnL0nGcM7MIq5YVEFlfvzuN60AFhGRD+ke9LOhpY9hfyhuZzdba9nQ0s+TG9pobh8gPTWZKxdVcMn8MgoyXbFu3gkpgEVE5KDhQIjmtg+6m4vjsLs55Di8ub2HJxva2ds7TEGmi5tX1nDBnBIyXIkTa4nTUhERGTeJ0N3sDYT5w5Yunmlsp8cToCo/nS9+ZAarZhSSkhx/d+gnogAWEZni4r27uW84wHMbO3ixuRNPIMzc8mw+d3YdS6rz4u6XhFOhABYRmaIOdDfv7R0mJz017rqb2/q8PN3YzqvvdxN2LCvqCrhiUQUzS7Ji3bSoUACLiEwxB7qbN7UPkpxkKM1xx9Wd5Pudgzy5oY11u/eTkmw4b3YJly8spyw3vn5BGCsFsIjIFHKgu9njD8XVZhqOtby7Zz9PbWhnS+cgmWnJfHxpJRfPLyM3Pb63uRwtBbCIyBRwZHdzvGymEQw7vL51H081ttHW56Moy8Vnz5zGebNLcKcmx7p540oBLCIyicVrd7PHH+KlTZ08u7GDvuEgtYUZfPn8mZwxvZDkONkqcrwpgEVEJql4nN3cM+Tn2aYOXt7chTcYZmFlLl86r4IFFTlx8YvBRFIAi4hMMvE4u3lv7zBPNrTxxrYeLJYzpxdy+aIK6ooyY920mFEAi4hMEkduphHr7mZrLZvaB3iyoZ31e/tIS0niovmlXLagLC5+KYg1BbCIyCQQ6W7ez7A/HPPuZsex/HlXL082tLG920OOO4VPLKvionmlZLsn54zm0VAAi4gksCNLBcbyzjIQcvjj+1083dhO54Cfshw3nzu7jnNnFeNKif34c7xRAIuIJKB42rt50BfkheZOfr+xgwFfiBnFmXzqo9Oon5ZP0hSZ0TwaCmARkQRiraV7yE/D3j48gXBMN9PoHvTxdGMHr2zpwh9yOL0mjysWVTCnLHvKzWgeDQWwiEiC8PhDNLcP0LrfS447lZLs2Iyn7tzn4cmGNt7e0YMxhrNnFnH5wnKqCzJi0p5EpQAWEYlzobDDzn0eNrUPkJKcFJPuZmstja39PNnQTlNrP+mpyVy2sJyPLSinINM1oW2ZLBTAIiJxylpL14CPhtZ+vIHYzG4OO5a3dvTwZEMbu3uGyc9I5VMrarhwbgkZLkXIWOjsiYjEoSF/iI2t/bT1+chNT6V4grub/aEwf9zSzdON7XQN+qnIc/OX505n1cwiUpM1ozkaFMAiInEkONLdvLl9gNTkJEpzJra7ecgX4vnmDp7b2MGgL8Sskiw+c8Y0Tp+WT5ImVkWVAlhEJA5Ya+kc8NHQ0o8vGKYgM21CixL0DPl5urGdlzdHZjQvrc7jqsUVzNaM5nGjABYRibFBX5CNbQO093vJS3dN6G5Rh+7RDHDWjEKuXFwxpWY0ewNhPP4gDlCaM3EbmSiARURiJBh22N41xJbOQVzJSZTlpE/Ye2/uGODJDe28u2f/wT2aL19YTlFW2oS1IVYcaxn2h/EFQ1ggNz2V+ZW5FGWnkZ02cbGoABYRmWDWWjr6fWxo7SMQciicoO5mx1re29PHkxva2NI5SFZaCtedXsUl8yf/Hs1hxzLkD+EPhTFASY6buRXZFGamke5KjkmbFMAiIhNowBdkY2s/Hf0+8jJc5LrHfw1tKOzwp+09PLmhjdY+L0VZLj57Zi3nzS7GnRqb8JkIwbDDoC9E2LEkGajMT6cyL528DFdc7E2tABYRmQCBkMO2riG2dQ3iSkmmLHf8u5t9wTAvb44UR+j1BKguyOCvzp/JGdMLYlotaTz5gmGG/CEca0lLSaauKIPSHDd5Ga4JndR2MhTAIiLjyFpLW5+XxtZ+gmFLYWbauBcoGPAGeW5jB883d+Dxh5lbns0XzqljcVXepJvRbK3FGwzj8YexWLLdqcwty6Y4201Oekpc/3sVwCIi46TfG6SppZ+uIT/5Gankpo9vd2/ngI+nG9t5ZUsXobClvjafqxZXMLMke1zfd6I51uLxh/AFwwAUZLqYWZxFUXYamRM4iWqsEqelIiIJwh8Kj3Q3D+FOSaZsnJe27OrxsGZDG2/t6CHJGM6ZWcQViyuozJu4WdXjLexYBn1BAmGHJBNZLlSZl0NhVlrCjmOPKYCNMV8DPg9YoBG41Vrri0bDREQSjeNYWvu8NLb24ThQNI7dzdZamtsHWLOhjYaWSHGEyydZcYQPJlE5JCcZKvPTqciNn0lUYzXqADbGVAJ3APOstV5jzCPAjcADUWqbiEjC6BsO0NDST6/HT35G2rgFhONY/ry7lyc3tLG920Nueio3LK/mormlCdX9eiz+UGQSVdiJ/0lUYzXW/60UIN0YEwQygLaxN0lEJHH4gmHe7xxkR/cQGa4USsdpM41g2OHVrd08taGdjgEfpTlpfO7sOs6dVZzwd4PeQJghfxALZKalcFppNqUJMIlqrEYdwNbaVmPMvcAewAs8b619PmotExGJY2HH0rJ/mKbWfgCKs93jUqxgOBDixeZOnm3qoM8bpK4ok69cOIsVtQXjPpt6vFhr8QTCeAMhAPIyXCysyqMoy0VW2uQO3UONpQs6H7gaqAP6gN8aYz5trf2vI153O3A7QE1NzehbKiISJ3o9ARr29tHvC5Kf4RqX8nz7hwM829jOi5u68AbDLKjM5a8WVzC/IichA8qxkZ2ofMHITlTFWW7mlGVTlBW7nahibSxd0B8FdlpruwGMMb8DzgIOC2Br7f3A/QD19fV2DO8nIhJT3kCYzR0D7O4ZJjMthZLs6M9ubu/38lRDO6++303YWlbWFXDlogqmF2dF/b3G26HbPyYbQ2lOGlUFueRnuBJ25nI0jSWA9wBnGGMyiHRBXwisjUqrRETiSNix7O7x0Nw2gDFQkh39Gr0793l4Yn0r7+zsJSXZcN7sYi5fWEFZ7sRV54mGUNhh0B8iGHZINoaqgsjM5fzM8ekpSGRjGQN+2xjzKPAuEALeY+ROV0Rksuge9LOhpQ+PP0RBhouUKIaItZZN7QM8sb6NhtbIUqIrF1fwsQVl5GUkzlKiYNhhwBvEsZbU5CRqCjIoy3WTPwlnLkfTmGZBW2u/DXw7Sm0REYkbHn+ITe0DtOz3ku2ObnezYy3v7tnPmvVtbO0aIic9lRuXV3PRvFIyXImxlOjQ5ULpqcnMKMmiLMdNbnpqwk4Om2iJ8T8tIjJBgmGHXfs8bGofICU5KardzWHH8sb2fazZ0EbLfi/FWWncuqqW804rSYilRJFCB0EcC5lpycwuzaZkCiwXGi8KYBERIt3BnQM+Glv78QXC5Ge6olYxKBByeOX9Lp7a0E73kJ+q/HS+dN4MzpxRGPdViYYDIYYDYRxryXGnMq8il+KRwvUK3bFRAIvIlDfgC9LcNkB7v5e8dBdF2dEpTj8cCPFCcyfPNHUw4A0yqySLz55Vy9KavHFZMxwNB6sLBUJgI2t0F1VmUZidRtYk2GkrnuhsisiUFanRO8i2riFSk5Moi9IuVn3DAZ5t6uCF5k68wTCLqnK5ekklc8uy4/Ku8ciNMQoyXcwqyaMoOy1hxqQTkc6siEw5jmNp749+jd7uQR9PNnxQDnBFXQFXL6mkrigzCq2OLmstHn8YbzASusVZbuaWRzbG0BrdiaEAFpEppW84QGNrPz1DAfKiVKN3b+8waza08cb2fRhjOHdWMVcuKqc8zsoBHqij6w2GMcZQkp3GvIrshC7pl8gUwCIyJRwomrCz20O6K5nSKNTo3do5yBMb2li3ez9pKUlcuqCcyxfGVznAI7eALM11s7Agl4JMF2kpCt1YUgCLyKQWdix7e4fZ2NaPtVCUnTamCVDWWhpa+lmzoY3m9gGy0lK47vQqLplfSrY7OpO3xupA6PqDYYyBshw31QV55GemKnTjiAJYRCatniE/DS399HuDFIxxK0THsbyzq5c1G9rYuc9DQaaLz5wxjQvmlMRF9+2hd7pJh4RuQebkKF4/GSmARWTSGQ5EdrHa0ztMdlrqmLqbQ2GH17bu48mGNtr7fZTnurn93OmcPbMo5nsbHy10awrzyM9Q6CYCBbCITBqhkV2smjsGSUkylGa7R73sxxcM8/LmLp5ubKfXE6CuKJOvXjiL5TGuw3u07mWFbmJSAItIwrPW0jXgo6G1H+8Yd7Ea8oV4bmMHv9/YwZA/xLzyHP7y3OksrMyN2Rrew2YvA+W56l6eDBTAIpLQBn1BNrYN0NHvI8edSvEod7HqGw7wTGM7L2zqxBd0WDYtn6sXVzCrNDvKLT45B9bpDgdDB2cvL1LoTioKYBFJSP5QmO1dQwd3sRrtOO++IT9PbmjjD1u6CDmWM6cXcvWSSmoKMqLc4hM7bEeqkXW68ysj63Q1e3nyUQCLSEJxHEtrn5emtn5CY9jFqr3fy5r1bby2dR8YOHdWEVcurqA8d2I3zziw9/KwP4QlsiPVvHJtjjEVKIBFJGHs90R2ser1+MnPSMOVfupdsXt6h3l8fStv7eghJcnw0XmlXLmonMKstHFo8bF5A5F6umApyHRxWmk+xdkK3alEASwicc8biOxitavHQ0ZqCqWjKJqwrWuIx9e3sm73ftypSVy5qIKPLSgjL2Pidq3yBcMM+YI4QG56KourcinOUcGDqUr/6yISt0Jhhz29wzS3DWAMFGWd2i5W1lo2tQ/w2Po2mlr7yUpL4fplVVwyv2zCSusFQg4DviCOtWSlpTC/MpeSHLdK+4kCWETij7WW7pFdrDz+EAUZLlJOYdMLay3r9/bx+PpW3u8cIi89lZtX1nDhnFLSXePfxRsKR0I35FjSU5OZU5ZNaa5bRezlMApgEYkrg74gm9oHaevzku1OoST75Gc3O9by5529PL6+lV09wxRlubh1VS3nnVYy7kt3wo5l0BckGHZITU6irjiTitx0ctNTFbpyVApgEYkLgZDD9u4h3u8cxJWcREl22kkHV8hxeGNbD09saKWtL7Jd5F+eO52zZxWNekOOk+FYy5AvhC8UJiXJUJWfQVV+OvkZrpjuliWJQQEsIjHlOJa2kWVFgVBkWVHySYZXMOzwx/e7WbO+je4hPzUFGdxxwSxW1o3fdpGHrdUFKvLSqSnIoCDz1LrJRRTAIhIz+z0Bmlr76RkOkJeeSm76yY3PHtin+amGNvYPB5lZksUtZ9WytCZv3Lp7fcEwg77IsqHCLBfzyvMpytYGGTJ6CmARmXAHlhXt3DdEhiuF0pMc5/X4Qzzf3MmzTe0M+iL7NH/pvJnMr8gZl+ANjkymCjuWHHcqi6pyKdGyIYkSXUUiMmFCYYfdPR42tQ+CgeJs90ktKxr0BXm2qYPnmjrwBsMsrc7jmqWVnDYO+zSHHRuZwRx2SEtN5rTSbMo0g1nGgQJYRMadtZbuQT8NrX14/OGTXlY04A3ydGM7zzd34A86rKgr4JqlldQWZka9fUMj1YZSkgzVBRlU5mkylYwvBbCIjKsBX5DmtgHa+3zkpqdSchLVivYPB3iqoZ2XNnUSCDmcOaOQa5ZUUh3lAgneQJihQBBspNrQ4sJItaFUTaaSCaAAFpFxcWS1orLcE4/z9gz5ebKhnZc3dxJ2LKtmFnHNkkoq8qJXIOHIcd0lVXmU5Li1B7NMOAWwiESV41ha93tpbOvDcTipakXdg37WbGjllS3dWAvnnlbE1UsqR11i8GhtGvAFCYQd0lKSmFmcRUV+Ojnu0dUOFokGBbCIRE3PkJ/Gln72ewMUZKSdcPepzgEfT6xv5dX3IyUBz59dzFWLKyg+hd2vjmc4EMLjD2EMVOZlUFOYQYHGdSVOKIBFZMyGAyE2tQ+wt9dLZloKZSeoVtTe5+Wx9a38ads+kqNcEjAYdhjwBglbS0GGi9nTImX+tF5X4o0CWERGLRh22LXPw6b2AVJOYvvIlv3DPPZeK2/u6CE1KYlLF5RzxaJy8sdYEtCxlgHvSBdzahKnlWVTnusmW13MEscUwCJyyqy1tPd7aWwdIBAMk5/pOu6ey7t7PDz2Xivv7OzFlZLEFQvLuXxRBbnpYwvISFH7IMYYKvPS1cUsCUUBLCKnpG84QFPrAPuG/OSmpx53ItOO7iEee6+Vtbv3k56azNVLKvnYwrIxTX4KOZEu5pBjyU1PZWl1PiW56mKWxKMAFpGT4guG2dIR2T4y3ZVy3BnKWzsH+d17razf20emK5nrTq/i0gVloy5Cb63F4w8zHAyRkmyoK8qkMj9jzHfQIrE0pgA2xuQBvwQWABa4zVr7ZhTaJSJxIhR22NM7THPbAOYE20e+3znI6nUtNLT2k5WWwg311Vw8v3TUeycHQg4D3gCWyPsuqs6lUFWHZJIY6x3wD4HnrLXXG2NcQHS3qRGRmLHW0jXgo7G1n+FAmPzjbB95aPDmuFO4aUUNF88rHdXmFo61DPpC+ENh3KnJzK3IoTw3ncxR3j2LxKtRX9HGmFzgXOAWAGttAAhEp1kiEksHto/s6PeR406l+BjbR27tHOTRd1toaIkE780ra/jo3NEFrz8UZsAbKfdXkZdObVGmJlTJpDaWXynrgG7g18aYxcA64CvWWk9UWiYiE84XDLOta4jtXUOkpSYfc5x3a+cgq99tYUNLP9nuFD61ooaLRnHHe+jyoUxXMgurcinP1baQMjWMJYBTgNOBv7bWvm2M+SFwF3D3oS8yxtwO3A5QU1MzhrcTkfESdix7e4fZ2NaPYy1FWUffPnJb1yCPrvsgeEfb1ewLhhn0BwFDVV460wozKMh0qdyfTCljCeAWoMVa+/bI548SCeDDWGvvB+4HqK+vt2N4PxGJMmst+4YCNLT0MegLHbMS0LauQVa/G5nVPNrgPexuNy2ZRZV5lOluV6awUQewtbbDGLPXGDPbWrsFuBBojl7TRGQ8DfqCbGofpK3PS7b76MuKtnUNsfrdFtbv7SMrLYWblldz8fyyUwrNQMih3xsELFX5GbrbFRkx1mmFfw3898gM6B3ArWNvkoiMpwNlArd2DeE6xvaRRwbvjcuruXheGemukwteOzKT2RcKk56qsV2RoxlTAFtr1wP10WmKiIynQ8sEhh0oOkqZwLEGbzDs0O8NYC2U5bqpK8qnMFMzmUWORgvrRKaAfSNlAvt9QfLTXR8qE7i9e4jV61p4byR4b1hezSWnELwefwhPIIQrOYnZZTlU5aePevMNkalC3yEik9iQP1ImsKV3mGx3KqVH1Nn9UPCews5V4ZEi98GwQ2GmiwWVhRRnp5Gsu12Rk6IAFpmEAiGH7d1DbO0cJDUlidIc92HjvLt7PPx2XQvrdu8nMy35lILXHwoz6A2CgdrCTGoKM7Uns8goKIBFJhHHsbT2eWlq6ycUthRkHn5H2rrfy2/X7eXtnb1kuJL5xLJIkYQTBa+1lkF/CF8wTKYrmUVVeZTluVWBSGQMFMAik0TPkJ/G1n76hgPkZ6ThSv9gnLej38fqd1v40/Z9pKUk8fGllVy2sPyE1YlCTmQJUdixlOa4Ob1Gk6pEokUBLJLgPCPjvHv3D5OVlkppTvrB57oH/Tz2Xgt/fL+blKQkLl9YzpWLK05Yj9cXDDPgC5KSZJheFOlmHm0pQRE5On1HiSSoQMhh574htnQMkpKcRGn2B+O8vZ4Aj69v5eXNXRjg4nllXL2kgrwM1zGPZ61lyB/CGwyTmZbM0uo8ynLTPzRjWkSiQwEskmAOjPNubOsneMQ4b783yJr1rbywqRPHgfNmF/PxpZUUZqUd93h93gChsKUkx81SdTOLTAgFsEgC6fUEaGztZ7/HT16Gi9z0yCSoIV+IJxva+P3GDgJhh3NnRYL3WNWMYKTYvS8IQG1hBtOKMk/YNS0i0aMAFkkAHn+ILZ0D7OnxkpmWcnCcdzgQ4pnGdp5p7MAXDHPmjEKuO72Kirz0Yx5rOBBiyB8iLTWJ+RU5VOSla4tIkRhQAIvEsUDIYVePh83tA6Qcsm+zLxjmuY0dPNXQhscfZkVtAdcvq6K6IOOoxzmwN7M3GCYvI5UVtQWU5Li1aYZIDCmAReKQ41ja+700tvYTCDsHx3kDIYcXmjtYs6GVAV+I02vyuH5ZNXVFmcc8Tp83QMixlOe6qS/OVyUikTihABaJM72eAE2t/fQOB8hLTyU33UXIcXhxUxe/e7eF/cNBFlbm8ollVcwqzT7qMSJFET4Y360tyiRb47sicUUBLBInhgMhNneMjPO6UijNduNYy5+27ePRdS10DPiYXZrNly+YxbzynKMe48D63dTkJOaWZVNVkKHxXZE4pQAWibFg2GHnPg+bOwZJSTKUZEeWDL27ez8Pr93Lnt5hphVkcOcls1lanXfU7mOPPzKxKsudzOnV+ZTluUlN1vpdkXimABaJEWsj47wNrf0Egg75mS5SkpLY1D7Aw3/ey5bOQcpy3Pz1BTM5Y3ohSUcE76ETqwozXSyqKqQo68M1fkUkPimARWJgvydAU1s/PUMBctNTyXW72LnPw8N/3sOGln7yM1L5/Nl1fGR2MSlJh9/JOtbSNxwkFHYoy3WzvLSA/IxUTawSSTAKYJEJ5A2E2dIxyK4eDxmuZEpz3LT1eXlkbaRCUVZaCjevrOHieWUf2gIyPDKjOexYagszqSvWxhkiiUwBLDIBgmGH3fs8bOoYJMlAcXYa+z0B7n91O398v5vU5CSuXVrJ5YvKP1QaMBR26BuZ0TyjOJNphZlkqjCCSMLTd7HIOLLW0tHvo7G1H18wTH6mi2F/mP96azcvNHcCcMn8Mq5eUvmhovaBkMOAN0BSkmFOWTbVmtEsMqkogEXGSf9wkKbWfrqH/OSmp5KcZHjsvVaeaWzHH3L4yKxirltWRdERhRIOLCVKS0liQVUelXmqSCQyGSmARaLMFwzzfucgO7s9uF3J5Ge4eKG5k8fXtzLkD7GyroBP1FdTecR+zd5AmAF/kExXZClReZ6bFC0lEpm0FMAiURIKO+zpHaa5bQBjoCDTxRvbe/jtur3sGwqwqDKXG5ZXM70467CvGw6EGPSFyHansKI2n9KcdO3RLDIFKIBFxshaS/egn4bWPjz+MPnpqWxsH+ChdyKbaNQVZfKX585gQWXuYV835A/h8QfJS3dx5oxCirWGV2RKUQCLjMGAL0hz2wAd/T5y3KkM+vz8/I87aG4foDQnjTsumMnKIzbRGPKFGA6GKMhwsbi6iOKsNK3hFZmCFMAio+ALhtnePcS2ziFcKUlYCz9/dTtv7+wlx53CrWfVcsGcksPGcAd9QYYDYYqyXCydVkShqhKJTGkKYJFTEHYsrfuHaWztxwKpyUn87r1WXt7URUqy4drTK7liYQXprshyIWstg/4Q3kCIkmw3y6apHKCIRCiARU7SviE/jS39I0uEknmhuYOnGtoJhh0umFPKdadXkpfhAg7dpzlESY6b5bXaLlJEDqcAFjkBjz/EpvYBWvYP405NZv2ePla/18qAN8jKugJuqK+mfGRJ0aEFEspy3awoLSA/0xXjf4GIxCMFsMgxHCwT2D5AcpJhR7eHR9ZG6vLOLc/mzotnM7MksqRIwSsip0oBLHKEw8oEhhw6+n3875/3smOfh+r8dL5xyWyWHFKXd8AXxBsIUZqTzoo6Ba+InBwFsMghDt0+ctAXZPW7razf20dhposvfmQG58wsOrhWd9AXZDgYpiQ7jeW1BRQoeEXkFCiARfhg+8gd3R4CYYfnmjr4w5YuMlKT+dSKGi6Z/0F5wCFfCE8gRHFWmmY1i8ioKYBlSgs7lr29w2xsi1QrenNHL09uaCPkWC6dX8bHl1aSPVJz1+MPMeQPUah1vCISBQpgmbL2DflpaOmjzxNkU8cAj65rYf9wkBV1Bdy0vIayXDcQ2at5yB8iL8PFquoiirIUvCIydmMOYGNMMrAWaLXWXjH2JomMr+FAiOa2yLKivb1eHl3Xwu7eYWYUZ3LHhbOYU5YDjFQn8gXJTU/lzOmFFGdry0gRiZ5o3AF/BdgE5EThWCLjJhh22LXPw+aOQToHfDy5oY0NLf0UZ0X2bD5jeiHGGPyhMH3eIFlpKaysK6A0x60iCSISdWMKYGNMFXA58D3gb6LSIpEos9bSOeCjsbWfzn4fL23u4pX3u0lPTebmlTVcPC8ywSoQcujz+klPTaZ+Wj7luSoLKCLjZ6x3wP8GfAPIHntTRKJv0BdkY9sAu3s8vLWjl2eb2gmGLJfMK+Pa0yMTrEJhh+5BH66UJJZW51GRl35YEQURkfEw6gA2xlwBdFlr1xljzjvO624HbgeoqakZ7duJnJJAyGF79xBb2gd4d08fTzW20+sJsKK2gBtXVFOem07Icdg35CfJwPyKHGoKM0lV8IrIBBnLHfAq4CpjzGWAG8gxxvyXtfbTh77IWns/cD9AfX29HcP7iZzQwV2sWvppbhvgifVtBydY/fX5M5lTnoPjWHqG/FjgtNIsaosySUtJjnXTRWSKGXUAW2u/CXwTYOQO+OtHhq/IROr3Bmlq6ae5Y4DnGttZt6ePoiwXXz5/JmfOKARg/3CAYNhhRnEWM4qzDpYNFBGZaFoHLAnPHwqzrWuIptZ+XtrUxcubu0hOMnyyvprLF5aTmmzo9wYJhBymFWYyqzSLzDRd+iISW1H5KWStfQV4JRrHEjlZjmNp7fPS0NLH61v38XRjOwO+EB85rZhP1ldTkOliyBdi/3CQ8tx05pTnkJueGutmi4gAugOWBNU/HKShtY+3dvSwZn0be/d7mV2azTcuncaM4iyGAyE6B7wUZqWxrLZEhRJEJO4ogCWh+ENh3u8Y5O2dvTzV0M76vZFx3jsumMUZ0wsIhB06B31kp6Vw1owi7V4lInFLASwJ4UB38zu7eni2oYNX3u8+bJw3yUT2dnalJrGsJp+KPG2iISLxTQEsca9/OMj6vft5prGDZ5sOH+fNS0+ld9hPkjFayysiCUUBLHHrwOzmZxvbeXx9Gy2HjPNOL8qkzxukx+NnZkkW04uzcKdqSZGIJA4FsMQday1tfV5e2tzF6nUtbGjppyjLxVcunMXKugIG/SG6h/xUF6RzWmkOWVpSJCIJSD+5JK4M+IL8eWcv//vOXv6wJbKe94b6ai5bWE7Icega9FOSk8bKugLyMjSzWUQSlwJY4kIw7LC1c5CH/7yXJza00Tcc5OyZRdy0ooastBT2ewORmc0zCynO0sxmEUl8CmCJqQOlAp9qaOehd/awvdtDbWEGX73wNGaUZLJ/OIA3GGJZdT4V+ZrZLCKThwJYYsbjD/Hm9n088MZu/rR9H5muFD5/dh0fOa2Yfl+QvuEgc8pyqC3MxJWimc0iMrkogGXChcKRUoG/fG0nzzS2MxwMc/G8Mq5fVkXYsfQOB6grymRWSbaKJYjIpKUAlgm1b8jPb9fu5b/e2k1rn4955Tl89qxairJcDHhDlOWmMa+iUHs2i8ikpwCWCeELhvnT1n389JVtvLunj8LMyLKipdV59HmDAKyapQlWIjJ1KIBlXDmOZUf3ED96aSvPNXdireXapZVctrAcbzDMcDAU2TpSE6xEZIpRAMu46fMEePCtXTz45m72DQWon5bPp8+owZWSzHAgxKzSbKYXZ5KWonFeEZl6FMASdYGQw6vvd3PfC+/T3D5ARZ6bb35sDtOLs/D4Q5RkpzG3PIdM7WAlIlOYfgJK1Fhr2d3j4V+ff5/fN3WQkmz49MppnDeniCFfCFdKEsumFas2r4gICmCJkuFAiAff3M39r+6g1xNg1YxCblxegzWWUNiyvLaA8tx0kjTOKyICKIBljMKO5e0dPfzTs5toah2gMs/N/7lsDpX5GYTCjjbSEBE5BgWwjFrngJd/fm4LT25oI8kYPrWihnNmFeEPOZTmpDGnTOO8IiLHop+OcsoCIYeH3tnDj1/eyr6hAGdML+CG+hqMsaS7kllRV0BhVlqsmykiEtcUwHLSrLVs2NvHd55q5r09fZTnurnr0jlU5qeTkmxYUJFLZZ7GeUVEToYCWE5K/3CQf/n9Zn67tgWAG+qrWDWzCICZJVnMKM7SOK+IyClQAMtxOY7lsfda+ZfnNtM56Gf5tHyur68iLSWZ0hw38ypyyHZr32YRkVOlAJZjer9zkLsfb+Ltnb2UZKfx9YtOo7owg8y0ZBZV5WnfZhGRMVAAy4f4AiHue2Erv3lzF461XH96JatmFZGalMTc8mymFWaSkqzuZhGRsVAAy2Fe2dLF3Y83sXe/l8VVedywvJoMVzK1hZnMLlN9XhGRaFEACwBdAz6+vWYjzzV1kJueypfPn8mMkkzyM1wsqsrT9pEiIlGmAJ7iHMfhP9/azfeff59Bf4iL55dx4dxiMl0pLKrMozJfy4pERMaDAngKa27r567fNdLQ0k9dUQZ/fcFMSrLdTC/OZFZpNu5UdTeLiIwXBfAU5A2E+Nffb+HBN3cfrFh0ek0uxTluFlXmkZuhZUUiIuNNATzFvLipk7sfb6K938eK2gKuWlxBYbZLu1iJiEwwBfAU0Tng5f881sSLm7ooynLx1xfMZGZJ1sE/aSnqbhYRmUgK4EnOcSwPvLGL77+wBV/Q4YpF5Zx7WjHluepuFhGJJQXwJLaxtZ87H22guX2AWSVZfKK+iqr8dBZW5qm7WUQkxkYdwMaYauBBoBSwwP3W2h9Gq2Eyer5gmP/vmU3851u7SU9N5i/OnMbS6lxmlmRrdrOISJwYyx1wCPhba+27xphsYJ0x5gVrbXOU2iaj8MqWLr75u0ba+32cNaOQyxaVU52fweKqXPIytJmGiEi8GHUAW2vbgfaRjweNMZuASkABHAO9Q37ufqKJpxs7KM5K444LZjK7NJsFVblU52eou1lEJM5EZQzYGFMLLAXejsbx5ORZa3l0XQv/z1PNDPlDXDK/lAvmlDC7NIc55epuFhGJV2MOYGNMFrAa+Kq1duAoz98O3A5QU1Mz1reTQ+zp8XDnow28vbOXaYUZfPG8GcwuzWZJdR6FWWmxbp6IiBzHmALYGJNKJHz/21r7u6O9xlp7P3A/QH19vR3L+0lEKOzw81d38KOXtmKB606v5JyZRSyoyqVWpQJFRBLCWGZBG+A/gE3W2vui1yQ5nsaWfv72t+t5v3OIeeU5XL+sknkVuSyozCUrTavKREQSxVh+Yq8CPgM0GmPWjzz2LWvtM2NulXyINxDmn5/bzINv7iIzLYXPnjmNM6YXsKQ6n7JcN5Hfh0REJFGMZRb064B+6k+AQ5cWrawr4JolFSyqzmNWSTauFHU3i4gkIvVZxrFeT4B/eKKJpxraKc5y8aXzZrCiroAl1Xla0ysikuAUwHHIWsvqd1v47lPNDPlCXDinhMsXlbO0Jp+aggyStaZXRCThKYDjTFuflzsf3cCftvVQU5DB7edMZ0VdIQsqc8hw6b9LRGSy0E/0OOE4lv98azf3PLuZsGO5anEFl84vZWmNJlmJiExGCuA4sHOfh799ZD3v7uljZkkWN9RXccb0Qk4ry1adXhGRSUoBHEOhsMP9r+7g317cSlISXL+sko8tKGNpTQEFmZpkJSIymSmAY2RzxwBfe3gDm9oHmFuWzU0razhreiF1xVmaZCUiMgUogCdYIOTww5fe52d/3IE7NYmblldz+aJyFlfnke1OjXXzRERkgiiAJ9D6vX38zcPr2bHPw+KqXG5aUcOqGUVU5qerXKCIyBSjAJ4A3kCYf/n9Zn7zxi6y3al89sxpXLm4gvkVuaS7NMlKRGQqUgCPsze393Dnoxto2e9l2bR8bl5Zw1kzCinN0dIiEZGpTAE8TgZ8Qb73VDMPr22hINPFF86p48rFFczW0iIREUEBPC5e3tzJ361uZN+Qn7NmFPKplTWcOb2Qwqy0WDdNRETihAI4ivZ7AvzDmiae3NBOaU4aXz5/JpctLGdmSRapyapaJCIiH1AAR8nvN3bwzd810jcc4PzZxXxqZQ0ragvJzdDSIhER+TAF8Bjt9wS4e6RkYHmum69ddBofW1BGXZE21BARkWNTAI/Bc00dfOuxyF3vBXNK+MwZNdTXFmhDDREROSEF8Cgcedf7NxedxhWLKqgpyNCGGiIiclIUwKfo0LveC+eWcMtZ0zi9poDMNJ1KERE5eUqNk9Q7ctf7dEM7FXlu7rx4NlcsrqAqP10baoiIyClTAJ+EA3e9/d4gH51bwm2rallaU6BtJEVEZNQUwMfR6wnwDyNjvZV5br5wzmlcs6SSslxtIykiImOjAD6G55ra+dZjTfR7g1w0r4TPrapjSU0+7lTd9YqIyNgpgI9w+F1vOn95znSuWVpBaW56rJsmIiKTiAL4EIfe9V4yr5Tbzq5lSU2+iieIiEjUKYA5fF1vZV46Xzx3OtecXklJtjvWTRMRkUlqygfwS5s6+bvVDewfDnLp/DJuXTWNpTUFuFJUPEFERMbPlA3gAV+Q//epZh5Z20JFrptvXDKb65ZVUaSSgSIiMgGmZAD/ads+7vztBtoHfHx0bgmfO7uOZdN01ysiIhNnSgXwcCDEPz+7md+8uZuS7DS+fvFsblherbteERGZcFMmgNft7uVvHtnA7p5hPnJaEV84ZzorpxeSmqy7XhERmXiTPoD9oTA/eGEr97+6nfwMF3/z0Vl86oxpuusVEZGYmtQB3NTaz9ceXs/WriFWzSjk8+fWcfbMYt31iohIzE3KAA6GHf79D9v58ctbyUxL4Y4LZvKZM2spztZdr4iIxIdJF8BbOwf520c20NDaz4raAj53dh0XzC3RXa+IiMSVSRPAYcfyq9d38q/Pb8GVnMSXPjKdW1bVUZKj3axERCT+jCmAjTGXAj8EkoFfWmvviUqrTtHuHg9f/+0G/rxrP0uq8/jCOXVcPL9Md70iIhK3Rh3Axphk4KfARUAL8GdjzBprbXO0Gnci1lr+6+09/NPTmwD4wjl13HJWLZX5GRPVBBERkVEZyx3wCmCbtXYHgDHmf4GrgQkJ4LY+L3+3uoHXtu5jQUUOt587nUsXlGs3KxERSQhjCeBKYO8hn7cAK498kTHmduB2gJqamjG83eH+4YmNvLOzl1vPquWWs2qZVpQZtWOLiIiMt3GfhGWtvR+4H6C+vt5G67j/eNU8ugf9zC3PwZ2qer0iIpJYxhLArUD1IZ9XjTw2IaryM6jSWK+IiCSosQyY/hmYZYypM8a4gBuBNdFploiIyOQ26jtga23IGPNl4PdEliH9ylq7MWotExERmcTGNAZsrX0GeCZKbREREZkytGZHREQkBhTAIiIiMaAAFhERiQEFsIiISAwogEVERGJAASwiIhIDCmAREZEYMNZGbXvmE7+ZMd3A7igesgjYF8XjTVU6j9Gh8xgdOo9jp3MYHdE4j9OstcVHe2JCAzjajDFrrbX1sW5HotN5jA6dx+jQeRw7ncPoGO/zqC5oERGRGFAAi4iIxECiB/D9sW7AJKHzGB06j9Gh8zh2OofRMa7nMaHHgEVERBJVot8Bi4iIJKSEDWBjzKXGmC3GmG3GmLti3Z5EYYzZZYxpNMasN8asHXmswBjzgjFm68jf+bFuZ7wxxvzKGNNljGk65LGjnjcT8aORa7PBGHN67FoeX45xHv/RGNM6ck2uN8Zcdshz3xw5j1uMMZfEptXxxxhTbYz5gzGm2Riz0RjzlZHHdU2epOOcwwm7HhMygI0xycBPgY8B84CbjDHzYtuqhHK+tXbJIdPr7wJestbOAl4a+VwO9wBw6RGPHeu8fQyYNfLnduD/TlAbE8EDfPg8Avxg5JpcMlJnnJHv6RuB+SNf8+8j3/sCIeBvrbXzgDOAvxo5X7omT96xziFM0PWYkAEMrAC2WWt3WGsDwP8CV8e4TYnsauA3Ix//Brgmdk2JT9baV4HeIx4+1nm7GnjQRrwF5BljyiekoXHuGOfxWK4G/tda67fW7gS2Efnen/Kste3W2ndHPh4ENgGV6Jo8acc5h8cS9esxUQO4Eth7yOctHP/EyQcs8LwxZp0x5vaRx0qtte0jH3cApbFpWsI51nnT9XnqvjzSNfqrQ4ZAdB5PgjGmFlgKvI2uyVE54hzCBF2PiRrAMnpnW2tPJ9Il9VfGmHMPfdJGpsVravwp0nkbk/8LzACWAO3A92PamgRijMkCVgNftdYOHPqcrsmTc5RzOGHXY6IGcCtQfcjnVSOPyQlYa1tH/u4CHiPShdJ5oDtq5O+u2LUwoRzrvOn6PAXW2k5rbdha6wC/4INuPZ3H4zDGpBIJjv+21v5u5GFdk6fgaOdwIq/HRA3gPwOzjDF1xhgXkYHxNTFuU9wzxmQaY7IPfAxcDDQROXefHXnZZ4EnYtPChHOs87YG+IuRmadnAP2HdAvKEY4Yi/w4kWsSIufxRmNMmjGmjsgEoncmun3xyBhjgP8ANllr7zvkKV2TJ+lY53Air8eUsXxxrFhrQ8aYLwO/B5KBX1lrN8a4WYmgFHgsct2RAvyPtfY5Y8yfgUeMMZ8jUq3qkzFsY1wyxjwEnAcUGWNagG8D93D08/YMcBmRSRrDwK0T3uA4dYzzeJ4xZgmR7tJdwF8CWGs3GmMeAZqJzFj9K2ttOAbNjkergM8AjcaY9SOPfQtdk6fiWOfwpom6HrUTloiISAwkahe0iIhIQlMAi4iIxIACWEREJAYUwCIiIjGgABYREYkBBbCIiEgMKIBFRERiQAEsIiISA/8/rgl3wqxTZb0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "model.plot_cumulative_hazard(figsize=(8,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
